Public Class Matrix3DSingle
    Implements IDisposable
    Public data(,,) As Single
    Public x, y, z As Integer
    Public Sub New()
        ReDim data(0, 0, 0)
        data(0, 0, 0) = 0
        x = 1
        y = 1
        z = 1
    End Sub
    Public Sub New(ByVal x As Integer, ByVal y As Integer, ByVal z As Integer)
        Me.x = x
        Me.y = y
        Me.z = z
        ReDim data(0 To x - 1, 0 To y - 1, 0 To z - 1)
        'automatically assigned 0
    End Sub
    Public Sub New(ByVal data(,,) As Single, ByVal i As Integer, ByVal j As Integer, ByVal k As Integer)
        Me.data = data
        Me.x = i
        Me.y = j
        Me.z = k
    End Sub
    Public Function Size() As Coordinate
        Return New Coordinate(Me.x, Me.y, Me.z)
    End Function
    Function Minimum() As Single
        Dim ans As Single = data(0, 0, 0)
        Dim i, j, k As Integer
        For i = 0 To x - 1
            For j = 0 To y - 1
                For k = 0 To z - 1
                    If data(i, j, k) < ans Then
                        ans = data(i, j, k)
                    End If
                Next
            Next
        Next
        Return ans
    End Function

    Function Maximum() As Single
        Dim ans As Single = data(0, 0, 0)
        Dim i, j, k As Integer
        For i = 0 To x - 1
            For j = 0 To y - 1
                For k = 0 To z - 1
                    If data(i, j, k) > ans Then
                        ans = data(i, j, k)
                    End If
                Next
            Next
        Next
        Return ans
    End Function

    Function TransverseSlice(ByVal zslice As Integer) As Single(,)
        Dim ans(Me.x, Me.y) As Single
        Dim i, j As Integer
        For i = 0 To Me.x - 1
            For j = 0 To Me.y - 1
                ans(i, j) = Me.data(i, j, zslice)
            Next
        Next
        Return ans
    End Function
    Function CoronalSlice(ByVal yslice As Integer) As Single(,)
        Dim ans(Me.x, Me.z) As Single
        Dim i, j As Integer
        For i = 0 To Me.x - 1
            For j = 0 To Me.z - 1
                ans(i, Me.z - 1 - j) = Me.data(i, yslice, j)
            Next
        Next
        Return ans
    End Function
    Function SagittalSlice(ByVal xslice As Integer) As Single(,)
        Dim ans(Me.y, Me.z) As Single
        Dim i, j As Integer
        For i = 0 To Me.y - 1
            For j = 0 To Me.z - 1
                ans(Me.y - 1 - i, Me.z - 1 - j) = Me.data(xslice, i, j)
            Next
        Next
        Return ans
    End Function
    Private Disposed As Boolean = False
    Private Overloads Sub Dispose(ByVal disposing As Boolean)
        ' Check to see if Dispose has already been called.
        If Not Me.disposed Then
            ' If disposing equals true, dispose all managed 
            ' and unmanaged resources.
            If disposing Then
                ' Dispose managed resources.
                x = Nothing
                y = Nothing
                z = Nothing
                data = Nothing
            End If

            ' Call the appropriate methods to clean up 
            ' unmanaged resources here.
            ' If disposing is false, 
            ' only the following code is executed.
        End If
        disposed = True
    End Sub
    ' Implement IDisposable.
    ' Do not make this method virtual.
    ' A derived class should not be able to override this method.
    Public Overloads Sub Dispose() Implements IDisposable.Dispose
        Dispose(True)
        ' This object will be cleaned up by the Dispose method.
        ' Therefore, you should call GC.SupressFinalize to
        ' take this object off the finalization queue 
        ' and prevent finalization code for this object
        ' from executing a second time.
        GC.SuppressFinalize(Me)
    End Sub

    Public Sub RotateAndTranslate(ByVal ThetaX As Single, ByVal ThetaY As Single, ByVal ThetaZ As Single, ByVal TranslateX As Single, ByVal TranslateY As Single, ByVal TranslateZ As Single)
        Dim Xcenter, Ycenter, Zcenter As Single
        Xcenter = Me.x / 2
        Ycenter = Me.y / 2
        Zcenter = Me.z / 2
        ThetaZ = -(ThetaZ * 3.14 / 180)
        ThetaX = -(ThetaX * 3.14 / 180)
        ThetaY = -(ThetaY * 3.14 / 180)

        Dim RotateMatrixX, RotateMatrixY, RotateMatrixZ, FinalRotateMatrix As New twoDMatrixLibrary.Matrix(4, 4)


        RotateMatrixX(0, 0) = 1
        RotateMatrixX(0, 1) = 0
        RotateMatrixX(0, 2) = 0
        RotateMatrixX(0, 3) = 0
        RotateMatrixX(1, 0) = 0
        RotateMatrixX(1, 1) = Math.Cos(ThetaX)
        RotateMatrixX(1, 3) = -1 * Math.Sin(ThetaX)
        RotateMatrixX(1, 3) = 0
        RotateMatrixX(2, 0) = 0
        RotateMatrixX(2, 1) = Math.Sin(ThetaX)
        RotateMatrixX(2, 2) = Math.Cos(ThetaX)
        RotateMatrixX(2, 3) = 0
        RotateMatrixX(3, 0) = 0
        RotateMatrixX(3, 1) = 0
        RotateMatrixX(3, 2) = 0
        RotateMatrixX(3, 3) = 1

        RotateMatrixY(0, 0) = Math.Cos(ThetaY)
        RotateMatrixY(0, 1) = 0
        RotateMatrixY(0, 2) = Math.Sin(ThetaY)
        RotateMatrixY(0, 3) = 0
        RotateMatrixY(1, 0) = 0
        RotateMatrixY(1, 1) = 1
        RotateMatrixY(1, 3) = 0
        RotateMatrixY(1, 3) = 0
        RotateMatrixY(2, 0) = -1 * Math.Sin(ThetaY)
        RotateMatrixY(2, 1) = 0
        RotateMatrixY(2, 2) = Math.Cos(ThetaY)
        RotateMatrixY(2, 3) = 0
        RotateMatrixY(3, 0) = 0
        RotateMatrixY(3, 1) = 0
        RotateMatrixY(3, 2) = 0
        RotateMatrixY(3, 3) = 1

        RotateMatrixZ(0, 0) = Math.Cos(ThetaZ)
        RotateMatrixZ(0, 1) = -1 * Math.Sin(ThetaZ)
        RotateMatrixZ(0, 2) = 0
        RotateMatrixZ(0, 3) = 0
        RotateMatrixZ(1, 0) = Math.Sin(ThetaZ)
        RotateMatrixZ(1, 1) = Math.Cos(ThetaZ)
        RotateMatrixZ(1, 3) = 0
        RotateMatrixZ(1, 3) = 0
        RotateMatrixZ(2, 0) = 0
        RotateMatrixZ(2, 1) = 0
        RotateMatrixZ(2, 2) = 1
        RotateMatrixZ(2, 3) = 0
        RotateMatrixZ(3, 0) = 0
        RotateMatrixZ(3, 1) = 0
        RotateMatrixZ(3, 2) = 0
        RotateMatrixZ(3, 3) = 1

        FinalRotateMatrix = RotateMatrixX * RotateMatrixY * RotateMatrixZ

        Dim TranslationmMatrix As New twoDMatrixLibrary.Matrix(1, 3)
        TranslationmMatrix(0, 0) = TranslateX
        TranslationmMatrix(0, 1) = TranslateY
        TranslationmMatrix(0, 2) = TranslateZ

        Dim temp(,,) As Single
        ReDim temp(x, y, z)
        Dim ii, jj, kk As Integer
        For k As Integer = 0 To z - 1
            For i As Integer = 0 To x - 1
                For j As Integer = 0 To y - 1
                    ii = CInt((FinalRotateMatrix(0, 0) * (i - Xcenter)) + (FinalRotateMatrix(0, 1) * (j - Ycenter)) + (FinalRotateMatrix(0, 2) * (k - Zcenter)) + Xcenter - TranslationmMatrix(0, 0))
                    jj = CInt((FinalRotateMatrix(1, 0) * (i - Xcenter)) + (FinalRotateMatrix(1, 1) * (j - Ycenter)) + (FinalRotateMatrix(1, 2) * (k - Zcenter)) + Ycenter - TranslationmMatrix(0, 1))
                    kk = CInt((FinalRotateMatrix(2, 0) * (i - Xcenter)) + (FinalRotateMatrix(2, 1) * (j - Ycenter)) + (FinalRotateMatrix(2, 2) * (k - Zcenter)) + Zcenter - TranslationmMatrix(0, 2))
                    If ii > 0 And jj > 0 And kk > 0 And ii <= x And jj <= y And kk <= z Then
                        temp(i, j, k) = data(ii - 1, jj - 1, kk - 1)
                    End If
                Next
            Next
        Next
        data = temp
    End Sub
    Public Function RotateAndTranslateFunction(ByVal ThetaX As Single, ByVal ThetaY As Single, ByVal ThetaZ As Single, ByVal TranslateX As Single, ByVal TranslateY As Single, ByVal TranslateZ As Single)
        Dim Xcenter, Ycenter, Zcenter As Single
        Xcenter = Me.x / 2
        Ycenter = Me.y / 2
        Zcenter = Me.z / 2
        ThetaZ = -(ThetaZ * 3.14 / 180)
        ThetaX = -(ThetaX * 3.14 / 180)
        ThetaY = -(ThetaY * 3.14 / 180)

        Dim RotateMatrixX, RotateMatrixY, RotateMatrixZ, FinalRotateMatrix As New twoDMatrixLibrary.Matrix(4, 4)


        RotateMatrixX(0, 0) = 1
        RotateMatrixX(0, 1) = 0
        RotateMatrixX(0, 2) = 0
        RotateMatrixX(0, 3) = 0
        RotateMatrixX(1, 0) = 0
        RotateMatrixX(1, 1) = Math.Cos(ThetaX)
        RotateMatrixX(1, 3) = -1 * Math.Sin(ThetaX)
        RotateMatrixX(1, 3) = 0
        RotateMatrixX(2, 0) = 0
        RotateMatrixX(2, 1) = Math.Sin(ThetaX)
        RotateMatrixX(2, 2) = Math.Cos(ThetaX)
        RotateMatrixX(2, 3) = 0
        RotateMatrixX(3, 0) = 0
        RotateMatrixX(3, 1) = 0
        RotateMatrixX(3, 2) = 0
        RotateMatrixX(3, 3) = 1

        RotateMatrixY(0, 0) = Math.Cos(ThetaY)
        RotateMatrixY(0, 1) = 0
        RotateMatrixY(0, 2) = Math.Sin(ThetaY)
        RotateMatrixY(0, 3) = 0
        RotateMatrixY(1, 0) = 0
        RotateMatrixY(1, 1) = 1
        RotateMatrixY(1, 3) = 0
        RotateMatrixY(1, 3) = 0
        RotateMatrixY(2, 0) = -1 * Math.Sin(ThetaY)
        RotateMatrixY(2, 1) = 0
        RotateMatrixY(2, 2) = Math.Cos(ThetaY)
        RotateMatrixY(2, 3) = 0
        RotateMatrixY(3, 0) = 0
        RotateMatrixY(3, 1) = 0
        RotateMatrixY(3, 2) = 0
        RotateMatrixY(3, 3) = 1

        RotateMatrixZ(0, 0) = Math.Cos(ThetaZ)
        RotateMatrixZ(0, 1) = -1 * Math.Sin(ThetaZ)
        RotateMatrixZ(0, 2) = 0
        RotateMatrixZ(0, 3) = 0
        RotateMatrixZ(1, 0) = Math.Sin(ThetaZ)
        RotateMatrixZ(1, 1) = Math.Cos(ThetaZ)
        RotateMatrixZ(1, 3) = 0
        RotateMatrixZ(1, 3) = 0
        RotateMatrixZ(2, 0) = 0
        RotateMatrixZ(2, 1) = 0
        RotateMatrixZ(2, 2) = 1
        RotateMatrixZ(2, 3) = 0
        RotateMatrixZ(3, 0) = 0
        RotateMatrixZ(3, 1) = 0
        RotateMatrixZ(3, 2) = 0
        RotateMatrixZ(3, 3) = 1

        FinalRotateMatrix = RotateMatrixX * RotateMatrixY * RotateMatrixZ

        Dim TranslationmMatrix As New twoDMatrixLibrary.Matrix(1, 3)
        TranslationmMatrix(0, 0) = TranslateX
        TranslationmMatrix(0, 1) = TranslateY
        TranslationmMatrix(0, 2) = TranslateZ

        Dim temp(,,) As Single
        ReDim temp(x, y, z)
        Dim ii, jj, kk As Integer
        For k As Integer = 0 To z - 1
            For i As Integer = 0 To x - 1
                For j As Integer = 0 To y - 1
                    ii = CInt((FinalRotateMatrix(0, 0) * (i - Xcenter)) + (FinalRotateMatrix(0, 1) * (j - Ycenter)) + (FinalRotateMatrix(0, 2) * (k - Zcenter)) + Xcenter - TranslationmMatrix(0, 0))
                    jj = CInt((FinalRotateMatrix(1, 0) * (i - Xcenter)) + (FinalRotateMatrix(1, 1) * (j - Ycenter)) + (FinalRotateMatrix(1, 2) * (k - Zcenter)) + Ycenter - TranslationmMatrix(0, 1))
                    kk = CInt((FinalRotateMatrix(2, 0) * (i - Xcenter)) + (FinalRotateMatrix(2, 1) * (j - Ycenter)) + (FinalRotateMatrix(2, 2) * (k - Zcenter)) + Zcenter - TranslationmMatrix(0, 2))
                    If ii > 0 And jj > 0 And kk > 0 And ii <= x And jj <= y And kk <= z Then
                        temp(i, j, k) = data(ii - 1, jj - 1, kk - 1)
                    End If
                Next
            Next
        Next
        Dim tempm As New Matrix3DSingle(temp, x, y, z)
        Return tempm
    End Function
    Default Public Property Item(ByVal x As Integer, ByVal y As Integer, ByVal z As Integer) As Single
        Get
            If (x < 0 Or y < 0 Or z < 0) Then Return Single.NaN
            If (x >= Me.x Or y >= Me.y Or z >= Me.z) Then Return Single.NaN
            Return data(x, y, z)
        End Get
        Set(ByVal value As Single)
            data(x, y, z) = value
        End Set
    End Property
    Default Public Property Item(ByVal xD As Double, ByVal yD As Double, ByVal zD As Double) As Single
        Get
            Return Me(Convert.ToInt16(Math.Round(xD)), Convert.ToInt16(Math.Round(yD)), Convert.ToInt16(Math.Round(zD)))
        End Get
        Set(ByVal value As Single)
            Me(Convert.ToInt16(Math.Round(xD)), Convert.ToInt16(Math.Round(yD)), Convert.ToInt16(Math.Round(zD))) = value
        End Set
    End Property
    Public Function ReScale(ByVal Data As Matrix3DSingle, _
    Optional ByVal MinimumValueOfRange As Single = 0, _
    Optional ByVal MaximumValueOfRange As Single = 255) As Matrix3DSingle

        Dim ReturnMatrix As Object
        ReturnMatrix = New Matrix3DSingle(Data.x, Data.y, Data.z)
        Dim Min As Single = Data.Minimum
        Dim Max As Single = Data.Maximum
        Dim Temp As Double
        For i As Integer = 0 To ReturnMatrix.x - 1
            For j As Integer = 0 To ReturnMatrix.y - 1
                For k As Integer = 0 To ReturnMatrix.z - 1
                    Temp = ((Data(i, j, k) - Min) / (Max - Min)) * (MaximumValueOfRange - MinimumValueOfRange) + MinimumValueOfRange
                        ReturnMatrix(i, j, k) = Temp

                Next
            Next
        Next
        Return ReturnMatrix
    End Function
    Public Function ReScaleToINT(ByVal Data As Matrix3DSingle, _
        Optional ByVal MinimumValueOfRange As Single = 0, _
        Optional ByVal MaximumValueOfRange As Single = 255) As Matrix3DInt16

        Dim ReturnMatrix As New Matrix3DInt16(Data.x, Data.y, Data.z)
        Dim Min As Single = Data.Minimum
        Dim Max As Single = Data.Maximum
        
        For i As Integer = 0 To ReturnMatrix.x - 1
            For j As Integer = 0 To ReturnMatrix.y - 1
                For k As Integer = 0 To ReturnMatrix.z - 1
                    If Data(i, j, k) <> 0 Then 'keep the noise same
                        ReturnMatrix(i, j, k) = Convert.ToInt16((Data(i, j, k) - Min) * ((MaximumValueOfRange - MinimumValueOfRange) / (Max - Min)) + MinimumValueOfRange)
                    End If
                Next
            Next
        Next

        Return ReturnMatrix
    End Function
    Sub InsertSlice(ByVal Begin0_End1 As Integer, ByVal C_T_S As String)
        Select Case C_T_S
            Case "C"
                Select Case Begin0_End1
                    Case 0
                        Dim temp As New Matrix3DSingle(Me.x, Me.y + 1, Me.z)
                        For a As Integer = 0 To Me.x - 1
                            For b As Integer = 0 To Me.y - 1
                                For c As Integer = 0 To Me.z - 1
                                    temp(a, b + 1, c) = data(a, b, c)
                                Next
                            Next
                        Next
                        data = temp.data
                    Case 1
                        Dim temp As New Matrix3DSingle(Me.x, Me.y + 1, Me.z)
                        For a As Integer = 0 To Me.x - 1
                            For b As Integer = 0 To Me.y - 1
                                For c As Integer = 0 To Me.z - 1
                                    temp(a, b, c) = data(a, b, c)
                                Next
                            Next
                        Next
                        data = temp.data
                End Select
            Case "T"
                Select Case Begin0_End1
                    Case 0
                        Dim temp As New Matrix3DSingle(Me.x, Me.y, Me.z + 1)
                        For a As Integer = 0 To Me.x - 1
                            For b As Integer = 0 To Me.y - 1
                                For c As Integer = 0 To Me.z - 1
                                    temp(a, b, c + 1) = data(a, b, c)
                                Next
                            Next
                        Next
                        data = temp.data
                    Case 1
                        Dim temp As New Matrix3DSingle(Me.x, Me.y, Me.z + 1)
                        For a As Integer = 0 To Me.x - 1
                            For b As Integer = 0 To Me.y - 1
                                For c As Integer = 0 To Me.z - 1
                                    temp(a, b, c) = data(a, b, c)
                                Next
                            Next
                        Next
                        data = temp.data
                End Select
            Case "S"
                Select Case Begin0_End1
                    Case 0
                        Dim temp As New Matrix3DSingle(Me.x + 1, Me.y, Me.z)
                        For a As Integer = 0 To Me.x - 1
                            For b As Integer = 0 To Me.y - 1
                                For c As Integer = 0 To Me.z - 1
                                    temp(a + 1, b, c) = data(a, b, c)
                                Next
                            Next
                        Next
                        data = temp.data
                    Case 1
                        Dim temp As New Matrix3DSingle(Me.x + 1, Me.y, Me.z)
                        For a As Integer = 0 To Me.x - 1
                            For b As Integer = 0 To Me.y - 1
                                For c As Integer = 0 To Me.z - 1
                                    temp(a, b, c) = data(a, b, c)
                                Next
                            Next
                        Next
                        data = temp.data
                End Select
        End Select
    End Sub
    Public Function Threshold(Optional ByVal MinimumValueOfThreshold As Single = 0, _
    Optional ByVal MaximumValueOfThreshold As Single = 100000) As Matrix3DSingle
        Dim ReturnMatrix As Object
        ReturnMatrix = New Matrix3DSingle(Me.data, Me.x, Me.y, Me.z)
        For i As Integer = 0 To ReturnMatrix.x - 1
            For j As Integer = 0 To ReturnMatrix.y - 1
                For k As Integer = 0 To ReturnMatrix.z - 1
                    If data(i, j, k) < MinimumValueOfThreshold Then ReturnMatrix(i, j, k) = 0
                    If data(i, j, k) > MaximumValueOfThreshold Then ReturnMatrix(i, j, k) = MaximumValueOfThreshold
                Next
            Next
        Next
        Return ReturnMatrix
    End Function
End Class